Brms_growth_width_Weibull

Go to the end of this file for the parameter estimates of the best fitted model on phenotype width

Trying to fit non-linear model with BRMS(Bayesian Regression Models using Stan)

widths <- read_csv("widths.csv")
widths

sample.list <- unique(widths$ID)

lapply(1:6, function(i) {
    ggplot(widths) +
      geom_line(aes(x=days, y=width, group=ID)) +
      facet_wrap_paginate(~ID, ncol = 6, nrow = 5,  page =i)
    }
)

Try weibull (see http://www.pisces-conservation.com/growthhelp/index.html?weibul.htm )

First set up the formula

weibull.bf1 <- bf(
  width ~ Hmax - (Hmax - Hmin) * exp(-(k*days)^delta), #Weibull model
  Hmax + Hmin + k + delta ~ 1, #general model, paramters do not vary for individuals
  nl=TRUE)

Priors. Hmin and Hmax use median at start and end dates.

stat.data <- function(data) {
  data %>% 
  group_by(days) %>%
  summarize(meadian=median(width),
            max=max(width),
            min=min(width),
            sd=sd(width))
}

stat.data(widths)

summary(fit1, waic=TRUE, R2=TRUE)

So with these priors it seems that the model can flip between having high and low Hmin and Hmax. This relates to the delta parameter flipping signs. So one possibility is to keep delta positive. Or more tightly constrain the priors on Hmax and Hmin.

summary(fit2, waic=TRUE, R2=TRUE)

 Family: gaussian(identity) 
Formula: width ~ Hmax - (Hmax - Hmin) * exp(-(k * days)^delta) 
         Hmax ~ 1
         Hmin ~ 1
         k ~ 1
         delta ~ 1
   Data: widths (Number of observations: 664) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1; 
         total post-warmup samples = 4000
    ICs: LOO = NA; WAIC = 5042.2; R2 = 0.32
 
Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept     59.78      0.71    58.43    61.17        344 1.00
Hmin_Intercept     29.69      5.14    16.08    36.98         50 1.10
k_Intercept         0.01      0.00     0.01     0.01         64 1.08
delta_Intercept     4.30      0.74     2.90     5.81         64 1.08

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma    10.78      0.29    10.21    11.34        170 1.05

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

plot(fit2)

pairs(fit2)

rescale k to get it more reasonable

weibull.bf2 <- bf(
  width ~ Hmax - (Hmax - Hmin) * exp(-((k/100)*days)^delta), #Weibull model
  Hmax + Hmin + k + delta ~ 1, #general model, paramters do not vary for individuals
  nl=TRUE)

summary(fit3, waic=TRUE, R2=TRUE)

 Family: gaussian(identity) 
Formula: width ~ Hmax - (Hmax - Hmin) * exp(-((k/100) * days)^delta) 
         Hmax ~ 1
         Hmin ~ 1
         k ~ 1
         delta ~ 1
   Data: widths (Number of observations: 664) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1; 
         total post-warmup samples = 4000
    ICs: LOO = NA; WAIC = 5040.41; R2 = 0.33
 
Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept     59.67      0.64    58.45    60.97       2415 1.00
Hmin_Intercept     31.18      5.50    17.27    38.14        692 1.00
k_Intercept         0.81      0.05     0.75     0.93        727 1.01
delta_Intercept     4.86      0.95     3.10     6.73        988 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma    10.74      0.30    10.17    11.35       2360 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

plot(fit3)

pairs(fit3)

What does the fit look like?

newdata <- data.frame(days=seq(min(widths$days), max(widths$days),1))
fit3.fitted <-  cbind(newdata, fitted(fit3, newdata)) %>% as.tibble() %>%
  rename(width=Estimate, lower.ci='2.5%ile', upper.ci='97.5%ile')

plot

pl <- ggplot(aes(x=days, y=width),data=NULL)
pl <- pl + geom_line(aes(group=ID), alpha=.2, data=widths)
pl + geom_line(color="skyblue", lwd=1.5, data=fit3.fitted)

Now try adding random effects for model parameters:

What parameters do we think might be interesting to allow to vary? Probably not Hmin. Try making a series of plots to see how varying delta or k affects things:

weibull.fn <- function(Hmax,Hmin,k,delta,days) {
    Hmax - (Hmax - Hmin) * exp(-((k/100)*days)^delta) #Weibull model
}
for(delta1 in seq(3,5,.25)) {
  tmp.widths <- weibull.fn(Hmax=62,
                           Hmin=20,
                           k=.7,
                           delta=delta1,
                           days=newdata$days)
  abc <- data.frame(newdata$days, tmp.widths)
  p <- ggplot(data=abc, aes(newdata$days, tmp.widths)) +
    geom_line() + ylim(0,100) + ggtitle("delta=",delta1)
  print(p)
}

for(k1 in seq(.5,1.5,.1)) {
  tmp.widths <- weibull.fn(Hmax=60,
                           Hmin=20,
                           k=k1,
                           delta=4,
                           days=newdata$days)
  abc <- data.frame(newdata$days, tmp.widths)
  p <- ggplot(data=abc, aes(newdata$days, tmp.widths)) +
    geom_line() + ylim(0,80) + ggtitle("k=",k1)
  print(p)
}

First try with only fixing Hmin

weibull.bf3 <- bf(
  width ~ Hmax - (Hmax - Hmin) * exp(-((k/100)*days)^delta), #Weibull model
  Hmax + k + delta ~ (1|ID), # vary for individuals
  Hmin ~ 1, # do not vary per individual
  nl=TRUE)

summary(fit5, waic=TRUE, R2=TRUE)

There were 12 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help.
See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

 Family: gaussian(identity) 
Formula: width ~ Hmax - (Hmax - Hmin) * exp(-((k/100) * days)^delta) 
         Hmax ~ (1 | ID)
         k ~ (1 | ID)
         delta ~ (1 | ID)
         Hmin ~ 1
   Data: widths (Number of observations: 664) 
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1; 
         total post-warmup samples = 10000
    ICs: LOO = NA; WAIC = 4943.66; R2 = 0.54
 
Group-Level Effects: 
~ID (Number of levels: 166) 
                    Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Hmax_Intercept)      3.94      1.07     1.34     5.81        526 1.01
sd(k_Intercept)         0.14      0.02     0.10     0.17       1611 1.00
sd(delta_Intercept)     0.64      0.59     0.02     2.22       4596 1.00

Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept     59.81      0.68    58.44    61.12       2333 1.00
k_Intercept         0.83      0.02     0.79     0.87       2558 1.00
delta_Intercept     7.54      0.51     6.52     8.52      10000 1.00
Hmin_Intercept     33.01      1.58    29.81    35.95       4256 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     8.84      0.33     8.23     9.51       2920 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

get fitted values

b <- widths %>% ungroup()
fit5.fitted <- cbind(b, fitted(fit5)) %>% as.tibble()
fit5.fitted

plot

plot.fitted <- function(fitted.data) {
   pl <- ggplot(fitted.data, aes(x=days)) +
     geom_line(aes(y=width),color="blue") +
     geom_line(aes(y=Estimate),color="red") +
     facet_wrap_paginate(~ID, ncol = 6, nrow = 5)
   pages <- n_pages(pl)
   lapply(seq_len(pages), 
          function(i) {
            ggplot(fitted.data, aes(x=days)) +
              geom_line(aes(y=width), color="blue") +
              geom_line(aes(y=Estimate), color="red") +
              facet_wrap_paginate(~ID, ncol = 6, nrow = 5, page=i)
          }
    )
}
plot.fitted(fit5.fitted)

plot fitted vs actual:

plot.fitted.actual <- function(fn) {
  
  r_squared <- summary(lm(Estimate ~ width, data=fn))$adj.r.squared
  r_squared <- round(r_squared, digits=4)

  fn %>%
    mutate(days=as.factor(days)) %>%
    ggplot(aes(x=width, y=Estimate, shape=days, color=days)) +
    geom_point() +
    geom_abline(intercept = 0, slope=1) +
    ggtitle(paste0("R2 = ",r_squared))

}

plot.fitted.actual(fit5.fitted)

total_sum_of_squared_residuals <- function(fn) {anova(lm(Estimate ~ width, data=fn))[2,2]}

cv

TSSR.fit <- total_sum_of_squared_residuals(fit5.fitted)
waic.fit <- round(waic(fit5)$waic, digits=2)
kfoldic.fit <- round(kfold(fit5)$kfoldic, digits=2)

widths.mean

widths.mean <- read_csv("widths.mean.csv")
widths.mean

lapply(1:6, function(i) {
       ggplot(aes(x=days, group=ID), data=widths.mean) +
         geom_line(aes(y=width), color="red") +
         geom_line(aes(y=width_raw), color="black") +
          facet_wrap_paginate(~ID, ncol = 6, nrow = 5,  page =i)
}
)

Try weibull (see http://www.pisces-conservation.com/growthhelp/index.html?weibul.htm )

stat.data(widths.mean) 

stat.data(widths) 

summary(fit3.mean, waic=TRUE, R2=TRUE)

 Family: gaussian(identity) 
Formula: width ~ Hmax - (Hmax - Hmin) * exp(-((k/100) * days)^delta) 
         Hmax ~ 1
         Hmin ~ 1
         k ~ 1
         delta ~ 1
   Data: widths.mean (Number of observations: 664) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1; 
         total post-warmup samples = 4000
    ICs: LOO = NA; WAIC = 4913.98; R2 = 0.48
 
Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept     65.46      0.84    63.98    67.26       2172 1.00
Hmin_Intercept     31.96      3.81    22.55    37.64       1081 1.00
k_Intercept         0.77      0.03     0.73     0.83       1355 1.00
delta_Intercept     4.60      0.77     3.20     6.16       1310 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     9.76      0.27     9.26    10.31       2722 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

plot(fit3.mean)

pairs(fit3.mean)

What does the fit look like?

newdata <- data.frame(days=seq(min(widths.mean$days), max(widths.mean$days),1))
fit3.mean.fitted <-  cbind(newdata, fitted(fit3.mean, newdata)) %>%
  as.tibble() %>%
  rename(width=Estimate, lower.ci='2.5%ile', upper.ci='97.5%ile')

plot

pl <- ggplot(aes(x=days, y=width),data=NULL)
pl <- pl + geom_line(aes(group=ID), alpha=.1, data=widths.mean)
pl + geom_line(color="skyblue", lwd=1.5, data=fit3.mean.fitted) +
  ggtitle("width.mean vs. days")

Now try adding random effects for model parameters:

What parameters do we think might be interesting to allow to vary? Probably not Hmin. Try making a series of plots to see how varying delta or k affects things:

weibull.fn <- function(Hmax,Hmin,k,delta,days) {
    Hmax - (Hmax - Hmin) * exp(-((k/100)*days)^delta) #Weibull model
}
for(delta1 in seq(3,8,.5)) {
  tmp.widths <- weibull.fn(Hmax=65,
                            Hmin=40,
                            k=1,
                            delta=delta1,
                            days=newdata$days)
   abc <- data.frame(newdata$days, tmp.widths)
  p <- ggplot(data=abc, aes(newdata$days, tmp.widths)) +
    geom_line() + ylim(0,80) + ggtitle("delta=",delta1)
  print(p)
}

for(k1 in seq(0,5,.5)) {
  tmp.width <- weibull.fn(Hmax=65,
                            Hmin=40,
                            k=k1,
                            delta=5,
                            days=newdata$days)
  abc <- data.frame(newdata$days, tmp.widths)
  p <- ggplot(data=abc, aes(newdata$days, tmp.widths)) +
    geom_line() + ylim(0,80) + ggtitle("k=",k1)
  print(p)
}

try with only fixing Hmin

weibull.bf3 <- bf(
  width ~ Hmax - (Hmax - Hmin) * exp(-((k/100)*days)^delta), #Weibull model
  Hmax + k + delta ~ (1|ID), # vary for individuals
  Hmin ~ 1, # do not vary per individual
  nl=TRUE)

summary(fit5.mean, waic=TRUE, R2=TRUE)

 Family: gaussian(identity) 
Formula: width ~ Hmax - (Hmax - Hmin) * exp(-((k/100) * days)^delta) 
         Hmax ~ (1 | ID)
         k ~ (1 | ID)
         delta ~ (1 | ID)
         Hmin ~ 1
   Data: widths.mean (Number of observations: 664) 
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1; 
         total post-warmup samples = 10000
    ICs: LOO = NA; WAIC = 3994.06; R2 = 0.92
 
Group-Level Effects: 
~ID (Number of levels: 166) 
                    Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Hmax_Intercept)      7.80      0.61     6.64     9.04       1874 1.00
sd(k_Intercept)         0.13      0.01     0.11     0.15       2589 1.00
sd(delta_Intercept)     2.49      0.36     1.85     3.25       1669 1.00

Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept     65.89      0.76    64.38    67.37       2669 1.00
k_Intercept         0.85      0.01     0.82     0.88       2056 1.00
delta_Intercept     5.43      0.42     4.63     6.28       4511 1.00
Hmin_Intercept     23.27      1.41    20.34    25.83       5660 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     3.78      0.22     3.38     4.23        993 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

get fitted values

b <- widths.mean %>% ungroup()
fit5.mean.fitted <- cbind(b, fitted(fit5.mean)) %>% as.tibble()
fit5.mean.fitted

plot

plot.fitted(fit5.mean.fitted)

plot fitted vs actual(modified):

plot.fitted.actual(fit5.mean.fitted)

plot fitted vs actual(raw):

plot.fitted.actual.raw <- function(fn) {
  
  r_squared <- summary(lm(Estimate ~ width_raw, data=fn))$adj.r.squared
  r_squared <- round(r_squared, digits=4)

  fn %>%
    mutate(days=as.factor(days)) %>%
    ggplot(aes(x=width_raw, y=Estimate, shape=days, color=days)) +
    geom_point() +
    geom_abline(intercept = 0, slope=1) +
    ggtitle(paste0("R2 = ",r_squared))

}

plot.fitted.actual.raw(fit5.mean.fitted)

TSSR on raw data set

total_sum_of_squared_residuals_raw <- function(fn) {anova(lm(Estimate ~ width_raw, data=fn))[2,2]}

TSSR.fit.mean.raw <- total_sum_of_squared_residuals_raw(fit5.mean.fitted)

cv

TSSR.fit.mean <- total_sum_of_squared_residuals(fit5.mean.fitted)
waic.fit.mean <- round(waic(fit5.mean)$waic, digits=2)
kfoldic.fit.mean <- round(kfold(fit5.mean)$kfoldic, digits=2)

widths.pmm

widths.pmm <- read_csv("widths.pmm.csv")
widths.pmm

 lapply(1:6, function(i) {
       ggplot(aes(x=days, group=ID), data=widths.pmm) +
         geom_line(aes(y=width), color="red") +   #red: modified data
         geom_line(aes(y=width_raw), color="black") +
          facet_wrap_paginate(~ID, ncol = 6, nrow = 5,  page =i)
}
)

Priors. Hmin and Hmax use median at start and end dates.

stat.data(widths)

stat.data(widths.mean)

stat.data(widths.pmm)

the stats of widths.mean and widths.pmm are very similar. so start the prior of width.mean.

weibull.bf3 <- bf(
  width ~ Hmax - (Hmax - Hmin) * exp(-((k/100)*days)^delta), #Weibull model
  Hmax + k + delta ~ (1|ID), # vary for individuals
  Hmin ~ 1, # do not vary per individual
  nl=TRUE)

summary(fit5.pmm, waic=TRUE, R2=TRUE)

 Family: gaussian(identity) 
Formula: width ~ Hmax - (Hmax - Hmin) * exp(-((k/100) * days)^delta) 
         Hmax ~ (1 | ID)
         k ~ (1 | ID)
         delta ~ (1 | ID)
         Hmin ~ 1
   Data: widths.pmm (Number of observations: 664) 
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1; 
         total post-warmup samples = 10000
    ICs: LOO = NA; WAIC = 4269.48; R2 = 0.88
 
Group-Level Effects: 
~ID (Number of levels: 166) 
                    Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Hmax_Intercept)      7.31      0.64     6.11     8.61       1578 1.00
sd(k_Intercept)         0.13      0.01     0.11     0.15       2101 1.00
sd(delta_Intercept)     2.43      0.38     1.74     3.21       1968 1.00

Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept     65.43      0.74    63.95    66.87       2204 1.00
k_Intercept         0.84      0.02     0.81     0.87       1922 1.00
delta_Intercept     5.84      0.46     4.97     6.78       4206 1.00
Hmin_Intercept     25.49      1.53    22.27    28.33       3241 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     4.76      0.25     4.29     5.27       1593 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

fitted values

a <- widths.pmm %>% ungroup()
fit5.pmm.fitted <- cbind(a, fitted(fit5.pmm)) %>% as.tibble()

plot fitted vs actual(modified)

plot.fitted.actual(fit5.pmm.fitted)

plot fitted vs actual(raw)

plot.fitted.actual.raw(fit5.pmm.fitted)

CV

0. TSSR on raw data set

TSSR.fit.pmm.raw <- total_sum_of_squared_residuals_raw(fit5.pmm.fitted)

1. TSSR

TSSR.fit.pmm <- total_sum_of_squared_residuals(fit5.pmm.fitted)

2. waic

waic.fit.pmm <- round(waic(fit5.pmm)$waic, digits=2)

3. kfoldic

#### which fitted model is the better?

create an empty list for models

0. TSSR on raw data set

cv.list[2,1] <- TSSR.fit.mean.raw
cv.list[3,1] <- TSSR.fit.pmm.raw

1. TSSR

cv.list[1,2] <- TSSR.fit
cv.list[2,2] <- TSSR.fit.mean
cv.list[3,2] <- TSSR.fit.pmm

2. CV

WAIC(Widely Applicable Information Criterion) an extension of the Akaike Information Criterion (AIC) that is more fully Bayesian.

K-fold CV Data are randomly partitioned into K subsets of equal size. Then the model is refit 10 times(default), each time leaving out one of the 10 subsets.

cv.list[1,3] <- waic.fit
cv.list[2,3] <- waic.fit.mean
cv.list[3,3] <- waic.fit.pmm

cv.list[1,4] <- kfoldic.fit
cv.list[2,4] <- kfoldic.fit.mean
cv.list[3,4] <- kfoldic.fit.pmm

cv.list

#### parameter estimates of best fitted model on phenotype width : take fit5.mean

read 
1) population (fixed) & group level (random) effects
2) parameters: Hmax, k & delta of each genotype

width.best.fitted <- fit5.mean$fit
dimnames(width.best.fitted)

$iterations
NULL

$chains
[1] "chain:1" "chain:2" "chain:3" "chain:4"

$parameters
  [1] "b_Hmax_Intercept"              "b_k_Intercept"                 "b_delta_Intercept"            
  [4] "b_Hmin_Intercept"              "sd_ID__Hmax_Intercept"         "sd_ID__k_Intercept"           
  [7] "sd_ID__delta_Intercept"        "sigma"                         "r_ID__Hmax[ID_1,Intercept]"   
 [10] "r_ID__Hmax[ID_10,Intercept]"   "r_ID__Hmax[ID_100,Intercept]"  "r_ID__Hmax[ID_101,Intercept]" 
 [13] "r_ID__Hmax[ID_102,Intercept]"  "r_ID__Hmax[ID_103,Intercept]"  "r_ID__Hmax[ID_104,Intercept]" 
 [16] "r_ID__Hmax[ID_105,Intercept]"  "r_ID__Hmax[ID_106,Intercept]"  "r_ID__Hmax[ID_107,Intercept]" 
 [19] "r_ID__Hmax[ID_109,Intercept]"  "r_ID__Hmax[ID_111,Intercept]"  "r_ID__Hmax[ID_112,Intercept]" 
 [22] "r_ID__Hmax[ID_113,Intercept]"  "r_ID__Hmax[ID_114,Intercept]"  "r_ID__Hmax[ID_116,Intercept]" 
 [25] "r_ID__Hmax[ID_117,Intercept]"  "r_ID__Hmax[ID_118,Intercept]"  "r_ID__Hmax[ID_119,Intercept]" 
 [28] "r_ID__Hmax[ID_12,Intercept]"   "r_ID__Hmax[ID_120,Intercept]"  "r_ID__Hmax[ID_121,Intercept]" 
 [31] "r_ID__Hmax[ID_122,Intercept]"  "r_ID__Hmax[ID_123,Intercept]"  "r_ID__Hmax[ID_124,Intercept]" 
 [34] "r_ID__Hmax[ID_125,Intercept]"  "r_ID__Hmax[ID_126,Intercept]"  "r_ID__Hmax[ID_127,Intercept]" 
 [37] "r_ID__Hmax[ID_128,Intercept]"  "r_ID__Hmax[ID_129,Intercept]"  "r_ID__Hmax[ID_13,Intercept]"  
 [40] "r_ID__Hmax[ID_130,Intercept]"  "r_ID__Hmax[ID_132,Intercept]"  "r_ID__Hmax[ID_133,Intercept]" 
 [43] "r_ID__Hmax[ID_134,Intercept]"  "r_ID__Hmax[ID_135,Intercept]"  "r_ID__Hmax[ID_136,Intercept]" 
 [46] "r_ID__Hmax[ID_137,Intercept]"  "r_ID__Hmax[ID_138,Intercept]"  "r_ID__Hmax[ID_139,Intercept]" 
 [49] "r_ID__Hmax[ID_14,Intercept]"   "r_ID__Hmax[ID_142,Intercept]"  "r_ID__Hmax[ID_143,Intercept]" 
 [52] "r_ID__Hmax[ID_144,Intercept]"  "r_ID__Hmax[ID_147,Intercept]"  "r_ID__Hmax[ID_148,Intercept]" 
 [55] "r_ID__Hmax[ID_15,Intercept]"   "r_ID__Hmax[ID_152,Intercept]"  "r_ID__Hmax[ID_153,Intercept]" 
 [58] "r_ID__Hmax[ID_154,Intercept]"  "r_ID__Hmax[ID_155,Intercept]"  "r_ID__Hmax[ID_156,Intercept]" 
 [61] "r_ID__Hmax[ID_157,Intercept]"  "r_ID__Hmax[ID_158,Intercept]"  "r_ID__Hmax[ID_16,Intercept]"  
 [64] "r_ID__Hmax[ID_160,Intercept]"  "r_ID__Hmax[ID_161,Intercept]"  "r_ID__Hmax[ID_163,Intercept]" 
 [67] "r_ID__Hmax[ID_165,Intercept]"  "r_ID__Hmax[ID_166,Intercept]"  "r_ID__Hmax[ID_169,Intercept]" 
 [70] "r_ID__Hmax[ID_17,Intercept]"   "r_ID__Hmax[ID_170,Intercept]"  "r_ID__Hmax[ID_171,Intercept]" 
 [73] "r_ID__Hmax[ID_172,Intercept]"  "r_ID__Hmax[ID_174,Intercept]"  "r_ID__Hmax[ID_175,Intercept]" 
 [76] "r_ID__Hmax[ID_176,Intercept]"  "r_ID__Hmax[ID_177,Intercept]"  "r_ID__Hmax[ID_18,Intercept]"  
 [79] "r_ID__Hmax[ID_180,Intercept]"  "r_ID__Hmax[ID_181,Intercept]"  "r_ID__Hmax[ID_186,Intercept]" 
 [82] "r_ID__Hmax[ID_189,Intercept]"  "r_ID__Hmax[ID_19,Intercept]"   "r_ID__Hmax[ID_191,Intercept]" 
 [85] "r_ID__Hmax[ID_192,Intercept]"  "r_ID__Hmax[ID_197,Intercept]"  "r_ID__Hmax[ID_198,Intercept]" 
 [88] "r_ID__Hmax[ID_199,Intercept]"  "r_ID__Hmax[ID_2,Intercept]"    "r_ID__Hmax[ID_20,Intercept]"  
 [91] "r_ID__Hmax[ID_201,Intercept]"  "r_ID__Hmax[ID_202,Intercept]"  "r_ID__Hmax[ID_207,Intercept]" 
 [94] "r_ID__Hmax[ID_208,Intercept]"  "r_ID__Hmax[ID_21,Intercept]"   "r_ID__Hmax[ID_212,Intercept]" 
 [97] "r_ID__Hmax[ID_214,Intercept]"  "r_ID__Hmax[ID_217,Intercept]"  "r_ID__Hmax[ID_223,Intercept]" 
[100] "r_ID__Hmax[ID_224,Intercept]"  "r_ID__Hmax[ID_227,Intercept]"  "r_ID__Hmax[ID_23,Intercept]"  
[103] "r_ID__Hmax[ID_231,Intercept]"  "r_ID__Hmax[ID_234,Intercept]"  "r_ID__Hmax[ID_24,Intercept]"  
[106] "r_ID__Hmax[ID_25,Intercept]"   "r_ID__Hmax[ID_26,Intercept]"   "r_ID__Hmax[ID_27,Intercept]"  
[109] "r_ID__Hmax[ID_28,Intercept]"   "r_ID__Hmax[ID_29,Intercept]"   "r_ID__Hmax[ID_3,Intercept]"   
[112] "r_ID__Hmax[ID_30,Intercept]"   "r_ID__Hmax[ID_32,Intercept]"   "r_ID__Hmax[ID_34,Intercept]"  
[115] "r_ID__Hmax[ID_35,Intercept]"   "r_ID__Hmax[ID_36,Intercept]"   "r_ID__Hmax[ID_37,Intercept]"  
[118] "r_ID__Hmax[ID_38,Intercept]"   "r_ID__Hmax[ID_39,Intercept]"   "r_ID__Hmax[ID_4,Intercept]"   
[121] "r_ID__Hmax[ID_40,Intercept]"   "r_ID__Hmax[ID_41,Intercept]"   "r_ID__Hmax[ID_42,Intercept]"  
[124] "r_ID__Hmax[ID_45,Intercept]"   "r_ID__Hmax[ID_47,Intercept]"   "r_ID__Hmax[ID_48,Intercept]"  
[127] "r_ID__Hmax[ID_49,Intercept]"   "r_ID__Hmax[ID_50,Intercept]"   "r_ID__Hmax[ID_51,Intercept]"  
[130] "r_ID__Hmax[ID_52,Intercept]"   "r_ID__Hmax[ID_53,Intercept]"   "r_ID__Hmax[ID_55,Intercept]"  
[133] "r_ID__Hmax[ID_56,Intercept]"   "r_ID__Hmax[ID_57,Intercept]"   "r_ID__Hmax[ID_58,Intercept]"  
[136] "r_ID__Hmax[ID_59,Intercept]"   "r_ID__Hmax[ID_6,Intercept]"    "r_ID__Hmax[ID_60,Intercept]"  
[139] "r_ID__Hmax[ID_61,Intercept]"   "r_ID__Hmax[ID_62,Intercept]"   "r_ID__Hmax[ID_63,Intercept]"  
[142] "r_ID__Hmax[ID_64,Intercept]"   "r_ID__Hmax[ID_65,Intercept]"   "r_ID__Hmax[ID_67,Intercept]"  
[145] "r_ID__Hmax[ID_68,Intercept]"   "r_ID__Hmax[ID_69,Intercept]"   "r_ID__Hmax[ID_7,Intercept]"   
[148] "r_ID__Hmax[ID_72,Intercept]"   "r_ID__Hmax[ID_73,Intercept]"   "r_ID__Hmax[ID_74,Intercept]"  
[151] "r_ID__Hmax[ID_75,Intercept]"   "r_ID__Hmax[ID_76,Intercept]"   "r_ID__Hmax[ID_77,Intercept]"  
[154] "r_ID__Hmax[ID_78,Intercept]"   "r_ID__Hmax[ID_79,Intercept]"   "r_ID__Hmax[ID_8,Intercept]"   
[157] "r_ID__Hmax[ID_80,Intercept]"   "r_ID__Hmax[ID_81,Intercept]"   "r_ID__Hmax[ID_82,Intercept]"  
[160] "r_ID__Hmax[ID_84,Intercept]"   "r_ID__Hmax[ID_85,Intercept]"   "r_ID__Hmax[ID_86,Intercept]"  
[163] "r_ID__Hmax[ID_87,Intercept]"   "r_ID__Hmax[ID_88,Intercept]"   "r_ID__Hmax[ID_89,Intercept]"  
[166] "r_ID__Hmax[ID_9,Intercept]"    "r_ID__Hmax[ID_90,Intercept]"   "r_ID__Hmax[ID_91,Intercept]"  
[169] "r_ID__Hmax[ID_92,Intercept]"   "r_ID__Hmax[ID_94,Intercept]"   "r_ID__Hmax[ID_95,Intercept]"  
[172] "r_ID__Hmax[ID_96,Intercept]"   "r_ID__Hmax[ID_98,Intercept]"   "r_ID__Hmax[ID_99,Intercept]"  
[175] "r_ID__k[ID_1,Intercept]"       "r_ID__k[ID_10,Intercept]"      "r_ID__k[ID_100,Intercept]"    
[178] "r_ID__k[ID_101,Intercept]"     "r_ID__k[ID_102,Intercept]"     "r_ID__k[ID_103,Intercept]"    
[181] "r_ID__k[ID_104,Intercept]"     "r_ID__k[ID_105,Intercept]"     "r_ID__k[ID_106,Intercept]"    
[184] "r_ID__k[ID_107,Intercept]"     "r_ID__k[ID_109,Intercept]"     "r_ID__k[ID_111,Intercept]"    
[187] "r_ID__k[ID_112,Intercept]"     "r_ID__k[ID_113,Intercept]"     "r_ID__k[ID_114,Intercept]"    
[190] "r_ID__k[ID_116,Intercept]"     "r_ID__k[ID_117,Intercept]"     "r_ID__k[ID_118,Intercept]"    
[193] "r_ID__k[ID_119,Intercept]"     "r_ID__k[ID_12,Intercept]"      "r_ID__k[ID_120,Intercept]"    
[196] "r_ID__k[ID_121,Intercept]"     "r_ID__k[ID_122,Intercept]"     "r_ID__k[ID_123,Intercept]"    
[199] "r_ID__k[ID_124,Intercept]"     "r_ID__k[ID_125,Intercept]"     "r_ID__k[ID_126,Intercept]"    
[202] "r_ID__k[ID_127,Intercept]"     "r_ID__k[ID_128,Intercept]"     "r_ID__k[ID_129,Intercept]"    
[205] "r_ID__k[ID_13,Intercept]"      "r_ID__k[ID_130,Intercept]"     "r_ID__k[ID_132,Intercept]"    
[208] "r_ID__k[ID_133,Intercept]"     "r_ID__k[ID_134,Intercept]"     "r_ID__k[ID_135,Intercept]"    
[211] "r_ID__k[ID_136,Intercept]"     "r_ID__k[ID_137,Intercept]"     "r_ID__k[ID_138,Intercept]"    
[214] "r_ID__k[ID_139,Intercept]"     "r_ID__k[ID_14,Intercept]"      "r_ID__k[ID_142,Intercept]"    
[217] "r_ID__k[ID_143,Intercept]"     "r_ID__k[ID_144,Intercept]"     "r_ID__k[ID_147,Intercept]"    
[220] "r_ID__k[ID_148,Intercept]"     "r_ID__k[ID_15,Intercept]"      "r_ID__k[ID_152,Intercept]"    
[223] "r_ID__k[ID_153,Intercept]"     "r_ID__k[ID_154,Intercept]"     "r_ID__k[ID_155,Intercept]"    
[226] "r_ID__k[ID_156,Intercept]"     "r_ID__k[ID_157,Intercept]"     "r_ID__k[ID_158,Intercept]"    
[229] "r_ID__k[ID_16,Intercept]"      "r_ID__k[ID_160,Intercept]"     "r_ID__k[ID_161,Intercept]"    
[232] "r_ID__k[ID_163,Intercept]"     "r_ID__k[ID_165,Intercept]"     "r_ID__k[ID_166,Intercept]"    
[235] "r_ID__k[ID_169,Intercept]"     "r_ID__k[ID_17,Intercept]"      "r_ID__k[ID_170,Intercept]"    
[238] "r_ID__k[ID_171,Intercept]"     "r_ID__k[ID_172,Intercept]"     "r_ID__k[ID_174,Intercept]"    
[241] "r_ID__k[ID_175,Intercept]"     "r_ID__k[ID_176,Intercept]"     "r_ID__k[ID_177,Intercept]"    
[244] "r_ID__k[ID_18,Intercept]"      "r_ID__k[ID_180,Intercept]"     "r_ID__k[ID_181,Intercept]"    
[247] "r_ID__k[ID_186,Intercept]"     "r_ID__k[ID_189,Intercept]"     "r_ID__k[ID_19,Intercept]"     
[250] "r_ID__k[ID_191,Intercept]"     "r_ID__k[ID_192,Intercept]"     "r_ID__k[ID_197,Intercept]"    
[253] "r_ID__k[ID_198,Intercept]"     "r_ID__k[ID_199,Intercept]"     "r_ID__k[ID_2,Intercept]"      
[256] "r_ID__k[ID_20,Intercept]"      "r_ID__k[ID_201,Intercept]"     "r_ID__k[ID_202,Intercept]"    
[259] "r_ID__k[ID_207,Intercept]"     "r_ID__k[ID_208,Intercept]"     "r_ID__k[ID_21,Intercept]"     
[262] "r_ID__k[ID_212,Intercept]"     "r_ID__k[ID_214,Intercept]"     "r_ID__k[ID_217,Intercept]"    
[265] "r_ID__k[ID_223,Intercept]"     "r_ID__k[ID_224,Intercept]"     "r_ID__k[ID_227,Intercept]"    
[268] "r_ID__k[ID_23,Intercept]"      "r_ID__k[ID_231,Intercept]"     "r_ID__k[ID_234,Intercept]"    
[271] "r_ID__k[ID_24,Intercept]"      "r_ID__k[ID_25,Intercept]"      "r_ID__k[ID_26,Intercept]"     
[274] "r_ID__k[ID_27,Intercept]"      "r_ID__k[ID_28,Intercept]"      "r_ID__k[ID_29,Intercept]"     
[277] "r_ID__k[ID_3,Intercept]"       "r_ID__k[ID_30,Intercept]"      "r_ID__k[ID_32,Intercept]"     
[280] "r_ID__k[ID_34,Intercept]"      "r_ID__k[ID_35,Intercept]"      "r_ID__k[ID_36,Intercept]"     
[283] "r_ID__k[ID_37,Intercept]"      "r_ID__k[ID_38,Intercept]"      "r_ID__k[ID_39,Intercept]"     
[286] "r_ID__k[ID_4,Intercept]"       "r_ID__k[ID_40,Intercept]"      "r_ID__k[ID_41,Intercept]"     
[289] "r_ID__k[ID_42,Intercept]"      "r_ID__k[ID_45,Intercept]"      "r_ID__k[ID_47,Intercept]"     
[292] "r_ID__k[ID_48,Intercept]"      "r_ID__k[ID_49,Intercept]"      "r_ID__k[ID_50,Intercept]"     
[295] "r_ID__k[ID_51,Intercept]"      "r_ID__k[ID_52,Intercept]"      "r_ID__k[ID_53,Intercept]"     
[298] "r_ID__k[ID_55,Intercept]"      "r_ID__k[ID_56,Intercept]"      "r_ID__k[ID_57,Intercept]"     
[301] "r_ID__k[ID_58,Intercept]"      "r_ID__k[ID_59,Intercept]"      "r_ID__k[ID_6,Intercept]"      
[304] "r_ID__k[ID_60,Intercept]"      "r_ID__k[ID_61,Intercept]"      "r_ID__k[ID_62,Intercept]"     
[307] "r_ID__k[ID_63,Intercept]"      "r_ID__k[ID_64,Intercept]"      "r_ID__k[ID_65,Intercept]"     
[310] "r_ID__k[ID_67,Intercept]"      "r_ID__k[ID_68,Intercept]"      "r_ID__k[ID_69,Intercept]"     
[313] "r_ID__k[ID_7,Intercept]"       "r_ID__k[ID_72,Intercept]"      "r_ID__k[ID_73,Intercept]"     
[316] "r_ID__k[ID_74,Intercept]"      "r_ID__k[ID_75,Intercept]"      "r_ID__k[ID_76,Intercept]"     
[319] "r_ID__k[ID_77,Intercept]"      "r_ID__k[ID_78,Intercept]"      "r_ID__k[ID_79,Intercept]"     
[322] "r_ID__k[ID_8,Intercept]"       "r_ID__k[ID_80,Intercept]"      "r_ID__k[ID_81,Intercept]"     
[325] "r_ID__k[ID_82,Intercept]"      "r_ID__k[ID_84,Intercept]"      "r_ID__k[ID_85,Intercept]"     
[328] "r_ID__k[ID_86,Intercept]"      "r_ID__k[ID_87,Intercept]"      "r_ID__k[ID_88,Intercept]"     
[331] "r_ID__k[ID_89,Intercept]"      "r_ID__k[ID_9,Intercept]"       "r_ID__k[ID_90,Intercept]"     
[334] "r_ID__k[ID_91,Intercept]"      "r_ID__k[ID_92,Intercept]"      "r_ID__k[ID_94,Intercept]"     
[337] "r_ID__k[ID_95,Intercept]"      "r_ID__k[ID_96,Intercept]"      "r_ID__k[ID_98,Intercept]"     
[340] "r_ID__k[ID_99,Intercept]"      "r_ID__delta[ID_1,Intercept]"   "r_ID__delta[ID_10,Intercept]" 
[343] "r_ID__delta[ID_100,Intercept]" "r_ID__delta[ID_101,Intercept]" "r_ID__delta[ID_102,Intercept]"
[346] "r_ID__delta[ID_103,Intercept]" "r_ID__delta[ID_104,Intercept]" "r_ID__delta[ID_105,Intercept]"
[349] "r_ID__delta[ID_106,Intercept]" "r_ID__delta[ID_107,Intercept]" "r_ID__delta[ID_109,Intercept]"
[352] "r_ID__delta[ID_111,Intercept]" "r_ID__delta[ID_112,Intercept]" "r_ID__delta[ID_113,Intercept]"
[355] "r_ID__delta[ID_114,Intercept]" "r_ID__delta[ID_116,Intercept]" "r_ID__delta[ID_117,Intercept]"
[358] "r_ID__delta[ID_118,Intercept]" "r_ID__delta[ID_119,Intercept]" "r_ID__delta[ID_12,Intercept]" 
[361] "r_ID__delta[ID_120,Intercept]" "r_ID__delta[ID_121,Intercept]" "r_ID__delta[ID_122,Intercept]"
[364] "r_ID__delta[ID_123,Intercept]" "r_ID__delta[ID_124,Intercept]" "r_ID__delta[ID_125,Intercept]"
[367] "r_ID__delta[ID_126,Intercept]" "r_ID__delta[ID_127,Intercept]" "r_ID__delta[ID_128,Intercept]"
[370] "r_ID__delta[ID_129,Intercept]" "r_ID__delta[ID_13,Intercept]"  "r_ID__delta[ID_130,Intercept]"
[373] "r_ID__delta[ID_132,Intercept]" "r_ID__delta[ID_133,Intercept]" "r_ID__delta[ID_134,Intercept]"
[376] "r_ID__delta[ID_135,Intercept]" "r_ID__delta[ID_136,Intercept]" "r_ID__delta[ID_137,Intercept]"
[379] "r_ID__delta[ID_138,Intercept]" "r_ID__delta[ID_139,Intercept]" "r_ID__delta[ID_14,Intercept]" 
[382] "r_ID__delta[ID_142,Intercept]" "r_ID__delta[ID_143,Intercept]" "r_ID__delta[ID_144,Intercept]"
[385] "r_ID__delta[ID_147,Intercept]" "r_ID__delta[ID_148,Intercept]" "r_ID__delta[ID_15,Intercept]" 
[388] "r_ID__delta[ID_152,Intercept]" "r_ID__delta[ID_153,Intercept]" "r_ID__delta[ID_154,Intercept]"
[391] "r_ID__delta[ID_155,Intercept]" "r_ID__delta[ID_156,Intercept]" "r_ID__delta[ID_157,Intercept]"
[394] "r_ID__delta[ID_158,Intercept]" "r_ID__delta[ID_16,Intercept]"  "r_ID__delta[ID_160,Intercept]"
[397] "r_ID__delta[ID_161,Intercept]" "r_ID__delta[ID_163,Intercept]" "r_ID__delta[ID_165,Intercept]"
[400] "r_ID__delta[ID_166,Intercept]" "r_ID__delta[ID_169,Intercept]" "r_ID__delta[ID_17,Intercept]" 
[403] "r_ID__delta[ID_170,Intercept]" "r_ID__delta[ID_171,Intercept]" "r_ID__delta[ID_172,Intercept]"
[406] "r_ID__delta[ID_174,Intercept]" "r_ID__delta[ID_175,Intercept]" "r_ID__delta[ID_176,Intercept]"
[409] "r_ID__delta[ID_177,Intercept]" "r_ID__delta[ID_18,Intercept]"  "r_ID__delta[ID_180,Intercept]"
[412] "r_ID__delta[ID_181,Intercept]" "r_ID__delta[ID_186,Intercept]" "r_ID__delta[ID_189,Intercept]"
[415] "r_ID__delta[ID_19,Intercept]"  "r_ID__delta[ID_191,Intercept]" "r_ID__delta[ID_192,Intercept]"
[418] "r_ID__delta[ID_197,Intercept]" "r_ID__delta[ID_198,Intercept]" "r_ID__delta[ID_199,Intercept]"
[421] "r_ID__delta[ID_2,Intercept]"   "r_ID__delta[ID_20,Intercept]"  "r_ID__delta[ID_201,Intercept]"
[424] "r_ID__delta[ID_202,Intercept]" "r_ID__delta[ID_207,Intercept]" "r_ID__delta[ID_208,Intercept]"
[427] "r_ID__delta[ID_21,Intercept]"  "r_ID__delta[ID_212,Intercept]" "r_ID__delta[ID_214,Intercept]"
[430] "r_ID__delta[ID_217,Intercept]" "r_ID__delta[ID_223,Intercept]" "r_ID__delta[ID_224,Intercept]"
[433] "r_ID__delta[ID_227,Intercept]" "r_ID__delta[ID_23,Intercept]"  "r_ID__delta[ID_231,Intercept]"
[436] "r_ID__delta[ID_234,Intercept]" "r_ID__delta[ID_24,Intercept]"  "r_ID__delta[ID_25,Intercept]" 
[439] "r_ID__delta[ID_26,Intercept]"  "r_ID__delta[ID_27,Intercept]"  "r_ID__delta[ID_28,Intercept]" 
[442] "r_ID__delta[ID_29,Intercept]"  "r_ID__delta[ID_3,Intercept]"   "r_ID__delta[ID_30,Intercept]" 
[445] "r_ID__delta[ID_32,Intercept]"  "r_ID__delta[ID_34,Intercept]"  "r_ID__delta[ID_35,Intercept]" 
[448] "r_ID__delta[ID_36,Intercept]"  "r_ID__delta[ID_37,Intercept]"  "r_ID__delta[ID_38,Intercept]" 
[451] "r_ID__delta[ID_39,Intercept]"  "r_ID__delta[ID_4,Intercept]"   "r_ID__delta[ID_40,Intercept]" 
[454] "r_ID__delta[ID_41,Intercept]"  "r_ID__delta[ID_42,Intercept]"  "r_ID__delta[ID_45,Intercept]" 
[457] "r_ID__delta[ID_47,Intercept]"  "r_ID__delta[ID_48,Intercept]"  "r_ID__delta[ID_49,Intercept]" 
[460] "r_ID__delta[ID_50,Intercept]"  "r_ID__delta[ID_51,Intercept]"  "r_ID__delta[ID_52,Intercept]" 
[463] "r_ID__delta[ID_53,Intercept]"  "r_ID__delta[ID_55,Intercept]"  "r_ID__delta[ID_56,Intercept]" 
[466] "r_ID__delta[ID_57,Intercept]"  "r_ID__delta[ID_58,Intercept]"  "r_ID__delta[ID_59,Intercept]" 
[469] "r_ID__delta[ID_6,Intercept]"   "r_ID__delta[ID_60,Intercept]"  "r_ID__delta[ID_61,Intercept]" 
[472] "r_ID__delta[ID_62,Intercept]"  "r_ID__delta[ID_63,Intercept]"  "r_ID__delta[ID_64,Intercept]" 
[475] "r_ID__delta[ID_65,Intercept]"  "r_ID__delta[ID_67,Intercept]"  "r_ID__delta[ID_68,Intercept]" 
[478] "r_ID__delta[ID_69,Intercept]"  "r_ID__delta[ID_7,Intercept]"   "r_ID__delta[ID_72,Intercept]" 
[481] "r_ID__delta[ID_73,Intercept]"  "r_ID__delta[ID_74,Intercept]"  "r_ID__delta[ID_75,Intercept]" 
[484] "r_ID__delta[ID_76,Intercept]"  "r_ID__delta[ID_77,Intercept]"  "r_ID__delta[ID_78,Intercept]" 
[487] "r_ID__delta[ID_79,Intercept]"  "r_ID__delta[ID_8,Intercept]"   "r_ID__delta[ID_80,Intercept]" 
[490] "r_ID__delta[ID_81,Intercept]"  "r_ID__delta[ID_82,Intercept]"  "r_ID__delta[ID_84,Intercept]" 
[493] "r_ID__delta[ID_85,Intercept]"  "r_ID__delta[ID_86,Intercept]"  "r_ID__delta[ID_87,Intercept]" 
[496] "r_ID__delta[ID_88,Intercept]"  "r_ID__delta[ID_89,Intercept]"  "r_ID__delta[ID_9,Intercept]"  
[499] "r_ID__delta[ID_90,Intercept]"  "r_ID__delta[ID_91,Intercept]"  "r_ID__delta[ID_92,Intercept]" 
[502] "r_ID__delta[ID_94,Intercept]"  "r_ID__delta[ID_95,Intercept]"  "r_ID__delta[ID_96,Intercept]" 
[505] "r_ID__delta[ID_98,Intercept]"  "r_ID__delta[ID_99,Intercept]"  "lp__"                         

r_ID__Hmax

datalist0 = list()
for (samples in sample.list) {
  ttt0 <- paste0("r_ID__Hmax[",samples,",Intercept]")
  dat0 <- extract(width.best.fitted, ttt0, permuted=TRUE)
  datalist0[[ttt0]] <- dat0
}
list.data0 = do.call(rbind, datalist0)
Hmax.ID <- data.frame(matrix(unlist(list.data0), nrow = 166, byrow = T))
row.names(Hmax.ID) <- paste0("r_ID__Hmax[",sample.list,",Intercept]")
head(Hmax.ID)

r_ID__k

datalist1 = list()
for (samples in sample.list) {
  ttt1 <- paste0("r_ID__k[",samples,",Intercept]")
  dat1 <- extract(width.best.fitted, ttt1, permuted=TRUE)
  datalist1[[ttt1]] <- dat1
}
list.data1 = do.call(rbind, datalist1)
k.ID <- data.frame(matrix(unlist(list.data1), nrow = 166, byrow = T))
row.names(k.ID) <- paste0("r_ID__k[",sample.list,",Intercept]")
head(k.ID)

r_ID__delta

datalist2 = list()
for (samples in sample.list) {
  ttt2 <- paste0("r_ID__delta[",samples,",Intercept]")
  dat2 <- extract(width.best.fitted, ttt2, permuted=TRUE)
  datalist2[[ttt2]] <- dat2
}
list.data2 = do.call(rbind, datalist2)
delta.ID <- data.frame(matrix(unlist(list.data2), nrow = 166, byrow = T))
row.names(delta.ID) <- paste0("r_ID__delta[",sample.list,",Intercept]")
head(delta.ID)

estimates of parameters

par.list = c("b_Hmax_Intercept", "b_k_Intercept", "b_delta_Intercept", "b_Hmin_Intercept", "sd_ID__Hmax_Intercept",    "sd_ID__k_Intercept","sd_ID__delta_Intercept","sigma", "lp__")
datalist3 = list()
for (pars in par.list) {
  dat3 <- extract(width.best.fitted, pars, permuted=TRUE)
  datalist3[[pars]] <- dat3
}
list.data3 = do.call(rbind, datalist3)
pars.est <- data.frame(matrix(unlist(list.data3), nrow = 9, byrow = T))
row.names(pars.est) <- par.list
pars.est

combine estimates of all parameters

pars.estimate.width.fitted <- rbind(pars.est, Hmax.ID, k.ID, delta.ID)
head(pars.estimate.width.fitted, 10)

write.csv(pars.estimate.width.fitted, file="pars.estimate.width.fitted.csv")